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We study two-level systems (2LS) coupled at different points to a one-dimensional waveguide in which one 
end is open and the other is either open (infinite waveguide) or closed by a mirror (semi-infinite). Upon injection 
of two photons (corresponding to weak coherent driving), the resonance fluorescence and photon correlations 
are shaped by the effective qubit transition frequencies and decay rates, which are substantially modified by 
interference effects. In contrast to the well-known result in an infinite waveguide, photons reflected by a single 
2LS coupled to a semi-infinite waveguide are initially bunched, a result that can be simply explained by stim¬ 
ulated emission. As the number of 2LS increases (up to 10 are considered here), rapid oscillations build up in 
the correlations that persist for a very long time. For instance, when the incoming photons are slightly detuned, 
the transmitted photons in the infinite waveguide are highly antibunched. On the other hand, upon resonant 
driving, incoherently reflected photons are mostly distributed within the photonic band gap and several sharp 
side peaks. These features can be explained by considering the poles of the single particle Green function in 
the Markovian regime combined with the time delay. Our calculation is not restricted to the Markovian regime, 
and we obtain several fully non-Markovian results. We show that a single 2LS in a semi-infinite waveguide can 
not be decoupled by placing it at the node of the photonic field, in contrast to recent results in the Markovian 
regime. Our results illustrate the complexities that ensue when several qubits are strongly coupled to a bus (the 
waveguide) as might happen in quantum information processing. 


I. INTRODUCTION 

The study of multiple photons confined in a one¬ 
dimensional (ID) waveguide interacting with local emitters 
(qubits) [ 1-5] has attracted a great deal of attention recently, 
and is now referred to by the term “waveguide quantum elec¬ 
trodynamics (QED)”. The ID geometry greatly limits the pos¬ 
sible propagating directions and hence increases the interfer¬ 
ence effects [4, 5] while decreasing the mode volume, which in 
turn enhances the coupling strength of qubits to the waveguide 
[6], The ID strong-coupling regime, where the light-matter 
interaction dominates over loss and dephasing, provides an 
excellent setting in which to investigate interesting quantum- 
optical effects theoretically [3-5, 7-34], observe such effects 
experimentally, [35-49], construct building blocks of quantum 
information processing and quantum computing [4, 7, 14, 50- 
58], and generate qubit-qubit entanglement [24, 59-64], 

A variety of artificial systems have been proposed and re¬ 
alized to implement light-matter interaction in ID, including 
superconducting qubits coupled to a microwave transmission 
line [35-49] or surface acoustic waves [65], and semiconduc¬ 
tor quantum dots coupled to either a metallic nanostmcture 
[66-68] or a photonic-crystal waveguide [69-71]. In addition 
to artificial atoms, waveguide QED can also be implemented 
using an ion trap [72], cold atoms trapped in [73] or near [74] 
an optical fiber, or single molecules doped in an organic crys¬ 
tal filled in a glass capillary [75]. In several of these systems, 
the coupling of the local emitter to the waveguide dominates 
by far all other emission or dephasing processes. 

The theoretical difficulty of waveguide QED lies in the 
fact that the waveguide photons are bidirectional while the 
qubits have arbitrary positions in the waveguide. A position- 
dependent phase factor is thus introduced even if the cou¬ 
pling strength for each qubit is the same. As a result, while 
a few photons scattering off one qubit (or multiple co-located 
qubits) has been extensively studied and exact solutions ex¬ 
ist [8-12, 14-17, 19, 21, 22, 29, 64, 76], for treating multi¬ 


ple qubits, the Markovian approximation has appeared nec¬ 
essary. Such Markovian multi-qubit, bidirectional waveguide 
calculations have been pursued recently using several theo¬ 
retical techniques: a Green function approach [77, 78], the 
master equation [60, 61, 79], input-output theory [26, 34], 
and the Lippmann-Schwinger (L-S) equation [24, 80], We 
note, however, one exception: an exact solution was ob¬ 
tained recently for two bidirectional photons scattering off 
two separated qubits [30], Furthermore, when entering the 
ultrastrong-coupling regime where the rotating-wave approx¬ 
imation (RWA) fails [81, 82], analytical treatments seem im¬ 
possible, and one has to use numerical methods such as matrix 
product states [32, 83, 84] to explore the many-body physics 
of photons. 

A single qubit in a semi-infinite waveguide is a more com¬ 
plex problem than for an infinite waveguide because of the 
delay in the reflection from the end and has therefore received 
considerable attention [49, 85-93]. Although an atom placed 
in front of a mirror in 3D open space has been studied both 
theoretically [85] and experimentally [94-97], the unconfined 
light in 3D makes the interference effect weak, and one there¬ 
fore expects a much stronger effect in ID. An exact solution 
for the wavefunction of the initially excited qubit can be de¬ 
rived by solving the delay-differential equation [85, 9-3], 
This solution demonstrates the complicated interference ef¬ 
fects caused by the mirror; if the distance to the mirror is large, 
non-Markovian effects come into play even for a single exci¬ 
tation (qubit or photon) [91-93]. Under the Markovian ap¬ 
proximation, the problem reduces to solving an ordinary dif¬ 
ferential equation [85, 93] which is far easier. The presence of 
the boundary in the semi-infinite case (i.e. the mirror) causes 
a modification of qubit frequencies and decay rates by modu¬ 
lating the stmcture of the photonic environment [49, 90], We 
are not aware of the existence of exact solutions for any cases 
of multi-photon scattering. 

In this paper, we consider N identical, equally-separated 
two-level systems (2LS) strongly coupled to an infinite or 
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semi-infinite waveguide (Fig. 1, the infinite waveguide has 
two open ends while the semi-infinite is closed by a perfect 
reflector on one end). Most of the results are obtained us¬ 
ing the Markovian approximation, which is checked by a full 
non-Markovian calculation in a few cases. It has been known 
since the introduction of the Dicke model [98] that interaction 
among the multiple 2LS can be induced through their cou¬ 
pling to bosonic modes, leading to sub- and super-readiance. 
In ID waveguides in particular, recent theoretical [26] and ex¬ 
perimental [46] studies of the power spectrum of two qubits 
coupled to an infinite waveguide clearly show that the qubit- 
qubit separation L modulates the effective resonant frequen¬ 
cies and decay rates, resulting in sub- and super-radiance. 
While it seems natural, then, to explore situations with many 
qubits, in fact discussion beyond two-qubit systems is limited 
in the literature [79, 80, 84, 99-101], Using the Lippmann- 
Schwinger equation, we show analytically that in the Marko¬ 
vian regime the collective behavior is encoded in the simple 
poles of the Green function. These poles reveal themselves in 
various measurable quantities such as the transmission spec¬ 
trum, time delay t, power spectrum S(u>) (resonance fluores¬ 
cence), and two-photon correlation functions ( 72 (f) (second- 
order coherence). The Markovian approximation reduces the 
number of poles from infinity to N [24] and so renders the 
problem tractable. Throughout the paper, we highlight a num¬ 
ber of common features of our results, such as rapidly oscil¬ 
lating two-photon correlations that persist for a long time, and 
the concentration of the reflected fluorescence within the pho¬ 
tonic band gap along with sharp side peaks. 

We point out an intriguing difference between the infinite 
and semi-infinite waveguides: while a single qubit coupled to 
the former can only reflect one photon at a time, giving rise to 
initial anti-bunching [7, 9, 102], when coupled to the latter it 
instead bunches the reflected photons. This can be explained 
simply by the stimulated emission. Another effect in a semi¬ 
infinite waveguide is the possibility of decoupling the waveg¬ 
uide from the 2LS by placing it at a node of the single-photon 
wavefunction when the qubit-mirror separation is small, as 
studied theoretically [90] and experimentally verified using 
superconducting qubits [49], We show that if the distance 
is large, however, the non-Markovian effects that come into 
play destroy this decoupling: our numerical non-Markovian 
calculation in the two excitation sector shows that the 2LS re¬ 
mains coupled to the waveguide because of oscillating non¬ 
trivial correlations. These two-photon features, to the best of 
our knowledge not addressed by previous ID studies [86-93] 
which mainly concern single-excitation properties, should be 
readily measurable using existing experimental technology. 

The rest of this paper is organized as follows: We first de¬ 
vote Sec. II to discussing the power spectrum of two pho¬ 
tons scattering off multiple distant qubits coupled to an infinite 
waveguide using the L-S equation. Since the power spectrum 
is a “first-order” quantity, one expects it to be easier to calcu¬ 
late and measure. In Sec. Ill we then move on to results for the 
second-order photon correlation . 92 (f)- To explain the long¬ 
time behavior of < 72 - we introduce the concept of time delay in 
Sec. IV. In Sec. V and VI, we turn to the discussion of power 
spectra and correlations for the semi-infinite waveguide. Some 
technical details are left for the appendices, including the de¬ 
tails of the L-S equation for both infinite and semi-infinite 



FIG. 1. (Color online) Schematic of the waveguide-QED system, in 
which equally separated, identical 2LS are coupled to a semi-infinite 
waveguide with one open end and another closed at x = 0. For an 
infinite waveguide with two open ends, the 2LS are instead placed 
symmetrically with respect to x = 0 to simplify the calculation. 


waveguides, the demonstration of the equivalence between the 
L-S equation and input-output theory at weak coherent driv¬ 
ing, and finally the two-photon transmission and reflection 
probabilities calculated using the L-S equation. 


H. MULTIPLE QUBITS IN AN INFINITE WAVEGUIDE: 
POWER SPECTRA 


Our starting point is the standard Hamiltonian used in 
waveguide QED [4, 7], consisting of a one-dimensional 
bosonic field coupled to discrete 2LS. After making the 
rotating-wave approximation (RWA) and extending the lim¬ 
its of the momentum integrals to infinity, one finds that the 
Hamiltonian in real space is (taking h = c = 1) 


H -Hqubit i I dx 


a t(x)-^ a R( x ) - al(x)-^a L (x) 


n y 

+ y F, V dx S(x- Xi) [a , a (x)ai- + a i+ a a {x)\ , 

A -1-T D J 


2—1 a—L,R _J' 


(i) 


where H qu b; t = ojq JV =1 cr i+ cr i-> a i± denotes the Pauli rais¬ 
ing (lowering) operator of the ;-th qubit with frequency ojq 
and position denotes the annihilation operator of right- 

(left-) going photons, and V is the coupling strength between 
the qubit and the photons. The decay rate for each qubit (to 
the waveguide) is T = 2V 2 . Throughout this paper we focus 
on the lossless limit, but loss could be simply introduced by 
tracing out an auxiliary waveguide or modifying the S-matrix 
elements [90, 103, 104], 

We calculate physical quantities by using the wavefunction 
\ip 2 ) = |^ 2 (fell ^ 2 ))rr of two incoming right-going photons 
(RR) with momenta k-\ and k- 2 - Throughout this work we fo¬ 
cus solely on identical incident photons: k\ = k 2 = E/2 
where E is the total input energy. The first physical quantity 
we consider is the power spectrum or resonance fluorescence, 
which is simply the Fourier transform of the first-order coher¬ 
ence. 


S a {uj) = J dte lult (ilj 2 \a'l{x 0 )a a (xo + t)\ip 2 ), (2) 
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where xq denotes the detector position (or equivalently, time) 
far away from the scattering region. S a (ui) is simply the spec¬ 
tral decomposition of the photons in the wavefunction \ip 2 ). 
Since the Hamiltonian H preserves the number of excitations 
(photon plus qubit), one can insert a one-particle identity op¬ 
erator, I dk\(f)i(k)) a ((f)-L(k)\, between the photon 

operators: 

S a (uj) = ]T / dk f dt e~ iuit 

a / =R,L* / 

x (i/j 2 \aa(xo)\(j)i(k)) a ' (4>i(k)\a a (xo + t)\ip 2 ), (3) 

where \<j>i(k)) a is the single-particle scattering eigenstate sat¬ 
isfying H\(j>\{k)) a = fc|0i(fc)) a with the incoming wave 
traveling in the a = L or R direction. The power spectrum 
follows by computing the matrix elements a '(^i(k)\a a {xo + 

t) 1 ^ 2 ). 

The two-photon wavefunction \ip 2 ) is obtained via the 
Lippmann-Schwinger equation following the procedure in 
Refs. [24, 80], The building blocks are the single-particle 
eigenstates \(f>i), the two-particle states \</> 2 ) formed from the 
direct product of two \<f>i), and the corresponding retarded 
Green function G R {E). In fact, \ip 2 ) can be written as [see 
Eq. (20) in Ref. 80] 

(^1 i k 2 )) ai ,Of2 — 102 (^-11 k 2 )) Qi ,a2 
N 

- E G R (E)\d l d l ) (G- 1 ).. (d J d J \Mh,k 2 )) ai , a2 . (4) 

i,j =1 


The second term contains all the nonlinearity and is often re¬ 
ferred to as the two-photon “bound state’’ [9, 12, 17], One 
can evaluate the desired matrix elements by inserting two- 
particle identity operators I 2 in the second term of Eq. (4) 
and performing the double momentum integral thereby intro¬ 
duced [see Eq. (23) in Ref. 80], This calculation is exact. For 
more than one qubit, making the Markovian approximation 
allows the integration to be done analytically. In this context, 
the Markovian approximation consists in replacing all factors 
of exp {ikL) that occur in the Green functions in Eq. (4) by 
exp(ifc 0 L), where k 0 = w 0 /c is the wavevector associated 
with wo and L = Xi +1 — Xi is the qubit-qubit separation 
[80]. In practice, we use a slightly modified expression for 
the power spectrum. 


S'r(w) =2Re^^ I dk 


dt t 


(5a) 


x {ip 2 \a'l(x 0 )\(l) 1 {k)) a '((l)i(k)\a R {x 0 +t)\4> 2 ), 


Sl(u) =2Re^] / dk 


dt / 


(5b) 


x (i/> 2 lal(x 0 )l</>i(k')) a '((/> 1 (k)la L (xo - t)\ijj 2 ). 


on resonance T=50% 



FIG. 2. (Color online) Normalized power spectra (resonance fluores¬ 
cence) of multiple qubits (from top to bottom: N = 1, 2, 3, 5,10) 
coupled to the infinite waveguide with koL = 7t/2 (separation 
L = Ao/4). For the first column, the incoming photons are on reso¬ 
nance, E/2 = loo = lOOT; for the second column the frequency is 
chosen such that the single photon transmission is 50% in each case. 
The total fluorescence (black solid line) is broken down into the re¬ 
flected (blue dashed) and transmitted (red dotted) components. The 
vertical lines indicate the real part of the poles. For N = 10 in the off- 
resonant case, the height of the central peak goes up to ~ 60, which 
is not shown for better visibility. The frequencies used in the second 
column are E/2T = {99.5, 99.29, 99.34, 99.43,99.48} (from top to 
bottom) [ 05]. 


where the former contains terms proportional to S(0)5(u) — 
E/2) because delta-normalized plane waves are used, and the 
latter is zero in the absence of the two-photon bound state and 
remains finite (for more discussion see Appendix D). Since the 
total incoherent/inelastic power spectrum is the sum of right- 
and left-going incoherent power spectra. 


The calculation of the matrix elements is given in Appendix 
A. 

After combining all pieces together, the resulting power 
spectrum can be divided into two parts, 

S a {u) = S“ herent M + 5“ mt ( w ), (6) 


^incoherent / 


t/, ,\ oincoherent/, ,\ i cmcoherent /, 

(W) — (w) + (wj, (/) 

one can normalize S'“ coherent (w) in terms of the incoherently 
scattered photon “flux” 


F 


incoherent 


duj 1 s ,inc ° herent 


(w). 


( 8 ) 
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FIG. 3. (Color online) Time delay t (top), single-photon transmission spectrum T = |t(fc )| 2 (middle), and poles of the transmission amplitude 
t(k ) (bottom) as a function of frequency. The system consists of 10 qubits coupled to an infinite waveguide with koL = n/A (left column) and 
koL = 7t/2 (right column). For the sake of clarity we show only the red-detuned side; for koL = n/2 the poles are symmetric with respect to 
the qubit frequency to o = 100T, while for koL = 7t/4 five poles are not shown. The vertical lines indicate the real parts of the poles. In panel 
(d) the black squares give the incident frequencies used in Fig. 8, the red dots give those used in Fig. 9, and the dashed line labels T = 50%. 
The dashed-dotted lines in panels (e) and (f) label the origin (T = 0). 


In this paper we omit the superscript “incoherent” for sim¬ 
plicity and focus on the inelastic power spectra normalized 
by F so that shapes and features can be readily compared. 
Furthermore, both on- and off-resonance cases are studied. 
In the off-resonant case, for a fair comparison of systems in 
which the number of qubits is different, we choose the inci¬ 
dent frequency such that (i) the single-photon transmission 
probability T is 50% in each case and (ii) it is the closest 
such frequency to the bare qubit frequency uj 0 (see discussion 
in Ref. [80]). Because we mostly focus on cases with small 
separation, k^L < 7r/2 (so L < Ao/4 with the wavelength 
Ao = 27 tc/wo). this choice leads to red-detuned incident fre¬ 
quencies, as will become clear in the following discussion. 

In interpreting the results, it will be useful to refer to the 
poles of the system, by which we mean the zeros of the denom¬ 
inator of the single-photon transmission or reflection ampli¬ 
tudes t(k) or r(k) [106]. Denote the poles by Zi = Ui — iTi/2 
with i = 1, 2, • • • , N (the factor of one half is in accordance 
with the definition of the decay rate 1'); then, the denomina¬ 
tor of t(k) and r{k) can be written as a polynomial of degree 
N, ( k — z±)(k — Z 2 ) ■ ■ ■ (k — 5/v). We will see that this in¬ 
deed gives us the effective qubit frequency and decay rate, as 
implied by the notation. In special cases the poles may be 
symmetrically arranged with respect to the u> = ujq line. This 
happens when k^L = n/2 because a wavefunction that is even 
about the middle of the interval between two adjacent qubits 
has the same magnitude at the site of those qubits as a wave- 
function that is odd. For other values of koL, the amplitude in 
these two cases is different, leading to an asymmetrical situa¬ 


tion in which there are superradiant and subradiant modes. 

The power spectra with qubit-qubit separation k 0 L = n/2 
are presented in Fig. 2. The result for the N = 1 case can be 
derived exactly and is given in Appendix A [Eq. (A9)]. In gen¬ 
eral, when the system is driven resonantly (E/2 = cuq) and the 
pole distribution is symmetric with respect to u > 0 (see Fig. 3 for 
a representative plot), both the transmitted and reflected power 
spectra are symmetric. In contrast, when the system is driven 
off-resonantly, neither the transmission nor the reflection fluo¬ 
rescence is symmetric. However, the total fluorescence (trans¬ 
mission + reflection) is still symmetric with respect to the in¬ 
cident frequency, indicating the conservation of energy and 
serving as a validity check on our calculation. In addition, 
thanks to the symmetric pole distribution for koL = n/2, the 
fluorescence when the incident photons are blue-detuned can 
be simply obtained by mirroring the red-detuned fluorescence 
with respect to ujq (data not shown). 

With regard to the dependence on the number of qubits, 
the main feature of the power spectra with resonant driving is 
that the photonic band gap develops, resulting in the decrease 
(increase) of transmission (reflection) fluorescence within the 
photonic band gap. In addition, many sharp side peaks appear 
around the photonic band gap, whose positions are roughly la¬ 
beled by the real parts of the poles { Cj, }. Since in general for 
large N the poles closer to uiq have smaller decay rates, we find 
that both the peak position and peak width could be explained 
by inspecting the poles’ real parts {Cji] and imaginary parts 
{Ti}, respectively. Finally, our two-qubit S(ui) agrees with 
the result obtained from input-output theory with weak coher- 
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FIG. 4. (Color online) Normalized power spectra (resonance fluo¬ 
rescence) of multiple qubits (from top to bottom: N = 2, 3, 5,10) 
coupled to the infinite waveguide with separation koL = 7t/4 (L = 
Ao/8). The incoming photon frequency is E/2 = ojo = 100T for 
the first column (on resonance) and is chosen such that T = 50% 
for the second column. The total fluorescence (black solid line) is 
broken down into the reflected (blue dashed) and transmitted (red 
dotted) components. The vertical lines indicate the real part of the 
poles. The frequencies used in the second column are E/2Y = 
{99.66, 99.73, 99.77, 99.78} (from top to bottom) [105]. 


ent driving [ 1 07], revealing the fact that two-photon scattering 
is the dominant process for weak driving. Further discussion 
is deferred to Appendix C. 

Furthermore, with slightly off-resonant driving the power 
spectra become sharply peaked. These sharp peaks reveal 
the existence of sub-radiant poles (with F < I ). Taking the 
N = 10 case as example [Fig. 2(j)], since the driving fre¬ 
quency is very close to the pole with the smallest F, that pole 
is highly excited and gives rise to the central peak with a very 
small width. The smaller peak on the right has the same F as 
the central peak and hence is visible too. Energy conservation 
then requires the smaller peak on the left to pop up as well. 
Thus, the fact that the poles largely determine the peak posi¬ 
tion and width is more transparent in the off-resonant cases, 
at least for those sub-radiant poles. We note that in any case 
transmission fluorescence is suppressed within the photonic 
band gap as expected. 

Next, we consider a smaller separation between the qubits, 
koL = 7t/ 4, see Fig. 4. For the resonant cases, the main dif¬ 
ference from the previous geometry koL = 7r/2 is that the 


large reflection fluorescence around u>o is reduced. Although 
it is still true that the reflection fluorescence is higher than the 
transmission fluorescence, the shape of the photonic band gap 
is different (red-detuned side is sharper than the blue-detuned 
side, connected to the asymmetric distribution of poles) mak¬ 
ing distinct peaks around {u>i} more visible. Energy conserva¬ 
tion implies, as before, that the total power spectrum is sym¬ 
metric with respect to u>o- The off-resonant sequence shows 
similar behavior to the kgL = it/2 case, with one minor dif¬ 
ference that the blue-detuned power spectra are different from 
the red-detuned spectra. In general the blue-detuned ones 
are much smoother because poles on the blue-detuned side 
(u>i > ojq) have larger decay rate F. Due to limited space 
we do not show them here. 


III. MULTIPLE QUBITS IN AN INLINITE WAVEGUIDE: 

PHOTON CORRELATIONS 

We now use the two-photon wavefunction |t/> 2 ) to calculate 
the second-order photon correlation function g 2 (t) (second or¬ 
der coherence), 

f f-j-\ = (&\ a t( x o)at,(xo + t)a Q (x 0 + tfaajxo)^) 
\(i/>2\ai(x 0 )a a (xo)\4>2}\ 2 

Since we are working in the two-photon sector, the numera¬ 
tor implies that g 2 is proportional to |(xo,xo + t|V > 2 )| 2 [80], 
We first show the cases with resonant driving and A: 0 L = tt/2 
or 7 t/4 in Fig. 5. Although the emergence of the two-photon 
bound state increases the probability for two photons to be 
transmitted [12], the L-S formalism in which an incoming 
plane-wave state is used gives an infinitesimally small correc¬ 
tion from the two-photon transmission (see Appendix D for 
details). Therefore, for simplicity we can ignore the transmis¬ 
sion 92 and focus on reflection g 2 for the resonant cases. 

It is well-known that a single 2LS cannot emit two photons 
at once because it can absorb only one photon at a time, so 
g 2 (0) = 0 in the reflection channel for N = 1 [ 7 , 9, 102], In 
contrast, we can see from Fig. 5 that adding more 2LS removes 
this limitation and allows g 2 ( 0) to be non-zero. The reason is 
that when one photon is trapped within the first 2LS, the other 
has a small chance to propagate to and be reflected by the next 
2LS, which in turn can cause the stimulated emission of the 
first photon. Thus, the probability of two photons coming out 
together is not fully suppressed, a scenario that is even more 
dramatic for the semi-infinite waveguide treated below. 

Secondly, note how oscillations build up and persist for a 
long time as N increases. We find that the frequency of long¬ 
time oscillations matches the difference between the incoming 
photon frequency E/2 = uq (the resonant frequency) and ui. t , 
the real part of the pole with the smallest decay rate Since 
the pole with the smallest decay rate occurs near the edge of 
the photonic band gap while the resonant frequency is near 
the middle of the gap, this low frequency scale should be oj ~ 
0.5F which is indeed what we observe. This makes sense since 
poles with larger decay rates have much less contribution to g 2 
at long time. In other words, we see the beating between the 
most sub-radiant pole and the driving frequency. 

We next discuss g 2 in the off-resonant cases. Because the 
N = 2 and 3 cases have been discussed in Ref. [80], here 
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FIG. 5. (Color online) Reflection g 2 of multiple qubits (from left to right: TV = 1, 2, 3, 5,10) coupled to an infinite waveguide with separation 
koL = 7t/2 (first row) and 7t/4 (second row). The two incoming photons are on resonance (E/2 = too = 100T). This comparison shows that 
the limitation <72 (0) = 0 for single 2LS is removed by adding more 2LS, and that a chain of few 2LS can cause long-time beating. 




FIG. 6. (Color online) In the off-resonant case, <72 of multiple qubits 
(left: N = 5; right: N = 10) coupled to an infinite waveguide with 
separation koL = 7r/2. The first row is for two transmitted photons 
and the second for two reflected ones. The gray, dashed line is the 
N = 1 result serving as a reference. The solid curve is calculated 
using the Markovian approximation while the full non-Markovian re¬ 
sult is given by the dots. The frequency of the incoming photons is 
chosen such that T = 50%, and the qubit frequency is u>o = 100T. 


FIG. 7. (Color online) In the off-resonant case, <72 of multiple qubits 
(left: N = 5; right: N = 10) coupled to an infinite waveguide with 
separation koL = 7r/4. The first row is for two transmitted photons 
and the second for two reflected ones. The gray, dashed line is the 
N = 1 result serving as a reference. The frequency of the incoming 
photons is chosen such that T=50%, and the qubit frequency is ojo = 
100T. 


we only present results for N = 5 and 10. For fc 0 L = 7t/2, 
they are shown in Fig. 6. It is known that, for k^L = 7t/2 
and N = 3, the transmission g> has a large initial bunching 
(g 2 > 1), while the reflection g 2 oscillates around the uncorre¬ 
lated value 1 [80]. It is striking that as N increases, the trans¬ 
mission correlations show antibunching (g 2 < 1) over a very 
long time, and the initial bunching is even diminished in the 
N = 10 case. The reflection g 2 continues to show a great 
deal of oscillation but in addition becomes highly bunched 
(g 2 > 1). The oscillation can be explained, as in the reso¬ 
nant case, by the beating between the most sub-radiant poles 
and the driving. 


We checked these results that use the Markovian approxi¬ 
mation against fully non-Markovian numerical results in a few 
cases. One of them is shown in Fig. 6(d). The agreement be¬ 
tween the two calculations (compare dots and solid line) is 
very good, showing that the Markovian approximation is rea¬ 
sonable for a qubit chain of moderate size. 

For k(, L = 7t/4 and off-resonant photons, the g 2 correla¬ 
tion is shown in Fig. 7. As in the N = 3 case [80], there is 
sharp initial bunching for both reflection and transmission. At 
non-zero t, the reflection g 2 shows bunching with irregular os¬ 
cillation while the transmission photons become strongly anti¬ 
bunched for a long time with little oscillation visible. The rea¬ 
son that g 2 of k$L = 7t/4 is very different from k^L = 7t/2 
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can be attributed to the highly asymmetric pole distribution. 
Take the TV = 10 case as an example for which the poles are 
shown in Fig. 3: there are two very close, sub-radiant poles 
that can contribute to the beating, and the beating frequency 
is small (one order of magnitude smaller than the tt/2 case) 
since we choose the incoming frequency to be red-detuned. 
The complicated interference effects result in highly nontrivial 
oscillations. We note that the oscillation is gradually washed 
out beyond Tt = 100 (data not shown). 

From the above results, g -2 is clearly very sensitive to the 
qubit-qubit separation L and the driving frequency (frequency 
of incoming photons). One may notice, however, that the res¬ 
onant cases with 7 t/ 2 and 7r/4 (Fig. 5) are somewhat more 
similar to each other and distinct from the off-resonant cases. 
Upon inspecting the polynomials giving rise to the poles for 
various system configurations, we find empirically that there 
is a general relation between the TV poles, 

1 iT 

o-y; (10) 

i =1 

that is, the average or “center of mass” of the poles coincides 
with the 2LS frequency and decay rate. This relation is inde¬ 
pendent of L and therefore provides a hand-waving explana¬ 
tion: upon resonant driving (E/2 = ujq), the incoming pho¬ 
ton frequency always matches the typical, average frequency 
of the excitations, leading to considerable absorption and ree¬ 
mission and and hence correlation. 

For the off-resonant case, we have chosen a particular value 
of the frequency for which T, the transmission, is 50%. There 
are, potentially, many such frequencies for a given system, and 
so we turn to comparing the behavior at these different points. 
As an example, we take the TV = 10, k^L = n/2 case. The 
chosen frequencies are labeled in Fig. 3(d); note that they are 
progressively further away from the resonance loq. The result 
is shown in Fig. 8. It is clear that the behavior is indeed some¬ 
what different for the five chosen frequencies. We first note 
that they all oscillate at roughly the same frequency due to the 
beating with the most sub-radiant poles. Secondly, the long¬ 
time structure of g 2 increases as the detuning becomes smaller, 
meaning that when driving very close to the frequency of the 
most sub-radiant pole [about 99.48r in this case; see Fig. 3(f)], 
a large time-scale sets in, leading to the long-time structure in 
<? 2 ■ As we shall see below, this is attributed to the large time 
delay associated with the most sub-radiant pole. 

In fact, if one drives very close to the most sub-radiant pole, 
the long-time structure is dramatic. We calculate three such 
frequencies [labeled in Fig. 3(d)] giving rise to T =20%, 50% 
(previously used), and 80%, respectively, and present the re¬ 
sult in Fig. 9. One can see that the long-time structure with 
off-resonant driving persists for more than Tt = 800 (much 
larger than the time of flight from one end of the array to the 
other without any obstacle, which is 97r/200r); in contrast, 
the time scale of the beating is almost invisible. Moreover, as 
one goes from T = 20% to T = 80% (approaching the sub¬ 
radiant pole), this long-time scale becomes larger, as if one 
“stretches” the g -2 curve. In the next section we employ the 
concept of time delay to explain this observation. 


IV. TIME DELAY 

The time delay (also known as the group delay) is a way 
to measure the time a wavepacket spends in passing through 
a scattering potential [108]. For a symmetric potential both 
transmitted and reflected wavepackets are characterized by a 
single time delay given by r(k) = dd^/dk in the general case 
and by 


t(/c) 


d 

dk 




e ik o L 


(ii) 


in the Markovian regime, where 0; c is the phase of the trans¬ 
mission amplitude t(k). 

A typical plot of the frequency dependence of the time de¬ 
lay is shown in Fig. 3 for TV = 10 with koL = 7 t/2 and 7r/4. 
It is clear that the position and width of the peaks in the time 
delay are precisely captured by, respectively, the real part {w*} 
and imaginary part {]%} of the poles. This means that as one 
approaches the sub-radiant poles, the time delay is greatly in¬ 
creased. In particular, for the most sub-radiant pole we find 
that the time delay t scales as TV 3 (fitting not shown), con¬ 
sistent with the finding by Tsoi & Law that the correspond¬ 
ing r scales as TV -3 [99]. Therefore, this feature explains the 
long-time structure of <72 discussed in the previous section: the 
“large time-scale” is contributed by the effective decay rate of 
the most sub-radiant pole. 

The flat structure of the time delay around uj 0 can also be 
explained. Within the photonic band gap, single photons are 
mostly reflected and hence spend much less time in the qubit 
array. Remarkably, we find empirically that the time delay at 
the resonant frequency is universal. 


T (fc=w 0 ) = ^, (12) 

independent of TV or L. This is consistent with the photon 
simply being reflected by the first qubit encountered. 

In short, the time delay is responsible for the long-time en¬ 
velope of <72 and it is directly connected to the simple poles 
of the system. In fact, for the off-resonant behavior in previ¬ 
ous sections, our choosing to work at the frequency closest to 
wo that satisfies T = 50% allowed us to take advantage of the 
associated long time-delay to examine nontrivial <y 2 behavior. 


V. MULTIPLE QUBITS IN A SEMI-INFINITE 
WAVEGUIDE: POWER SPECTRA 

We now turn to the case of a semi-infinite waveguide and 
study how the presence of a mirror (the boundary) changes 
the response of the system. As in the infinite waveguide case 
above, we first focus on the power spectra (fluorescence). The 
TV = 1 case has been analyzed by Koshino & Nakamura us¬ 
ing the Heisenberg-Langevin equation (equivalent to the input- 
output theory) at both weak and strong coherent driving [90], 
We find that our L-S approach gives the same result as theirs 
in the weak driving limit (see Appendix C), which hence vali¬ 
dates our calculation. 

Two changes in the calculation must be made for the semi- 
infinite case (see Appendix B for details). First, formally 
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FIG. 8. (Color online) of 10 qubits coupled to an infinite waveguide with separation koL = n/2. The first (second) row is for two transmitted 
(reflected) photons. The driving frequencies are chosen such that T = 50% and are labeled as black squares on the single-photon transmission 
spectrum in Fig. 3(d). The qubit frequency is u>o = 100T. 



FIG. 9. (Color online) Transmission g 2 of 10 qubits in an infinite 
waveguide with separation koL = tv /2. The driving frequencies are 
chosen such that the single particle transmission is (a) 20%, (b) 50%, 
and (c) 80% [labeled by red dots in Fig. 3(d)], close to the most sub¬ 
radiant pole. The qubit frequency is ojq = 100T. 


the Hamiltonian Eq. (1) remains the same, but the integration 
range is modified to be from negative infinity to zero. Accord¬ 
ingly, when solving for the single-particle eigenstate \vf>i(k)), 
a boundary condition f jv(fc) + r/v(fc) = 0 has to be imposed. 
We stress that in contrast to Koshino & Nakamura’s approach 
[90], here the boundary condition is imposed at the wavefunc- 
tion level rather than the Hamiltonian level, but the results 
agree exactly. Secondly, as there is only one incoming and 
outgoing channel, the summation over the incident direction 
a = {R, L} must be dropped. As a result, adding a mirror ac¬ 
tually reduces the number of matrix elements to be calculated. 


With the qubit-mirror separation defined to be |ccjv | = a, the 
Markovian approximation can be employed straightforwardly 
by replacing exp (ika) by exp(ifcoa), as done in the infinite 
waveguide case. 

In light of the discussion of the infinite waveguide case, we 
consider the case where the qubit-qubit separation is fixed at 
koL = tv 1% allowing the distribution of poles to be symmet¬ 
ric for certain values of a, and the qubit-mirror separation is 
varied to see how the mirror modifies the fluorescence. We 
focus mostly on the case of one and two qubits, commenting 
on the N = 10 results only at the end of this section. For one 
or two qubits and k^a = tv/2 or 7r/4, results are presented 
in Fig. 10 and Fig. 11. First, note that since the reflection 
fluorescence is the total fluorescence, the spectrum is always 
symmetric with respect to the incident frequency E/2. Sec¬ 
ond, since the single-photon reflection probability is always 
one, the way we chose the off-resonant driving frequency for 
the infinite case is no longer possible; instead, we have studied 
properties at fixed detunings. 

For N = 1 the results are similar to those of the infi¬ 
nite waveguide (cf. Figs. 2 and 4). Resonant driving gives 
a Lorentzian-like fluorescence, and off-resonant driving splits 
the Lorentzian peak into two. The condition for resonance is, 
of course, controlled by the single pole in this N = 1 case. 
The main difference here compared to the infinite waveguide 
case is that the pole is modulated by the qubit-mirror separa¬ 
tion a: 

r 

u! = u>o — — sin(2fc 0 a), T = T [1 — cos(2fc 0 a)]. (13) 

Thus, the effective frequency and decay rate of the qubit can 
be changed. These relations agree with those of Ref. [90] for 
a hard-wall boundary condition (9b = 7t/2 therein), and are 
responsible for the shift in the peak in the fc 0 a = 7r/4 case 
shown in Fig. 10(b). 

The spectrum changes dramatically compared with the infi¬ 
nite waveguide case when N > 2. For N = 2, the expressions 
for the poles are much more complicated than in the infinite 
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FIG. 10. (Color online) Normalized power spectra of one or two 
qubits coupled to a semi-infinite waveguide with qubit-qubit sepa¬ 
ration koL = 7r/2 (for N = 2). The qubit-mirror separation is 
koa = 7t/2 in (a) and (c), and n/4 in (b). The qubit frequency is 

ojo = ioor. 


waveguide case, and we do not reproduce them here. How¬ 
ever, for the special case k^L = 7r/2 we find that the poles can 
be simplified to 


5i, 2 (a) =w 0 -y±^v / 1- 2e 2ik ° a . (14) 

From this expression one can see that the “center of mass” is 
wo — iT/2 —it is not affected by the mirror. The two poles 
circulate this center in an elliptical trajectory on the complex 
plane as a changes, in contrast to the infinite waveguide case 
where the two poles circulate in a perfect circle as L changes 
[24, 99], 

For the case koa = 7r/2, the two poles have the same decay 
rate and the spectra are symmetric between red and blue detun¬ 
ing, so only the red-detuned case is shown in Fig. 10(c). When 
the driving frequency is far detuned, there are four peaks, with 
the inner two higher and the outer two lower, similar to that of 
the (total) power spectrum in the infinite case (not shown). The 
main difference is that here there are nodes (at which S = 0) 
between the inner and the outer peaks, one of which is fixed 
at the bare qubit frequency wo- As the driving frequency ap¬ 
proaches either of the poles, the two inner peaks merge into 
one (a process similar to that seen in Fig. 2). Next, when the 
driving frequency is tuned between the poles, both nodes start 
to be shifted and lifted, and do not touch down to zero again 
until the driving is on resonance (E/2 = wq). 



FIG. 11. (Color online) Normalized power spectra of two qubits 
coupled to a semi-infinite waveguide with qubit-qubit separation 
koL = 7r/2 and qubit-mirror separation koa = 7r/4. The driving 
frequencies for each plot are E/2Y = 96.5, 99,100,101, and 103.5 
(from bottom to top). The black ticks label the position of the two 
poles (and the qubit frequency too = 100T is at the center), and the 
blue arrow indicates the incident frequency. 


On the other hand, for the koa = 7r/4 case the decay rates of 
the two poles are different, resulting in a sharper (flatter) spec¬ 
trum on the red- (blue-) detuned side. To illustrate the drasti¬ 
cally varying structure of the fluorescence, we show in Fig. 11 
results for five incoming photon frequencies: substantially 
red-detuned, slightly red-detuned, likewise for blue-detuned, 
and finally on resonance. Starting from substantially red- 
detuned driving ( E/2 = 96.5T), the four peaks and two nodes 
are still visible, but the right node is red-shifted away from to o, 
presumably due to the asymmetric poles. As the frequency of 
the incoming photons is increased, the merging process hap¬ 
pens but with one difference from the koa = tt/2 case: the 
outer peaks disappear completely. For driving in between the 
poles, the main peak splits. In contrast to the /,: (l a = 7r/2 case, 
when the driving approaches the blue-detuned pole, instead of 
merging the two peaks actually shrink, and a single larger peak 
emerges between them. Finally, as the driving becomes sub¬ 
stantially blue-detuned, the larger peak again splits into two, 
with the outer peaks and the nodes appearing. 

We note a special case in this progression: atE/2 = 99.5T, 
the entire spectrum of inelastic scattering disappears and both 
photons are reflected elastically. The reason is that in the 
steady state the wavefunctions for the two qubits differ by a 
phase 7r. Together with the phases picked up during propaga¬ 
tion, it results in a precise destructive interference killing the 
photon-photon bound state. 

In Fig. 12(a) we show the power spectrum for a representa¬ 
tive N = 10 case. Compared to the infinite waveguide case 
(cf. Fig. 2), note the better defined photonic band gap behavior 
around the 2LS resonant frequency and the sharper modulation 
on the sides. This comes about because the mirror effectively 
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FIG. 12. (Color online) A representative case of (a) S(to) and (b) 
772 for 10 qubits coupled to the semi-infinite waveguide with koL = 
koa = 7t/2. The system is driven resonantly {E/2 = too = 100T). 
The vertical lines indicate the real parts of the poles. 


doubles the number of qubits that the photons see, leading to 
finer and stronger interference effects. 

In short, adding a mirror changes drastically the spectrum 
of inelastic scattering by two qubits and brings in another way 
to modulate the distribution of the poles. More generally, this 
will be the case for changing the boundary condition on the 
semi-infinite waveguide. For superconducting qubits coupled 
to a microwave transmission line, while physically moving 
the qubit in situ is normally not feasible, changing the bound¬ 
ary condition continuously with a magnetic field is readily ac¬ 
complished by adding a SQUID to the end of the waveguide 
[90, 109], 


VI. MULTIPLE QUBITS IN A SEMI-INFINITE 
WAVEGUIDE: PHOTON CORRELATIONS 

Finally, let us turn to results for photon correlations in a 
semi-infinite waveguide. We first concentrate on the single- 
2LS case. Because properties are controlled by a single pole 
[with frequency and decay rate give in Eq. (13)], g 2 will be the 
same for driving frequencies equally detuned (either blue- or 
red-) from to. The result is shown in Fig. 13 for both koa = 
7t/ 2 and 7r/4 (a = Ao/4 or Ao/8, respectively). 

A striking difference from the infinite waveguide case is that 
< 72 ( 0 ) is no longer zero; instead, it indicates bunching in all 
four cases shown. This can be explained by stimulated emis¬ 
sion: since the first photon is captured by the 2LS, the second 
photon passes through to the wall and is reflected back. Be¬ 
cause of the short distance (time-of-flight = 2 a/c~ 7 r/wo < C 
1 /F), when the second photon revisits the 2LS, the first photon 



FIG. 13. (Color online) <72 of a single qubit coupled to a semi-infinite 
waveguide. The qubit-mirror separation is (a),(c) koa = 7 r /2 and 
(b),(d) 7 t/ 4. The frequency of the incoming photons is (a),(b) reso¬ 
nant with the 2LS ( E/2 = too) and (c),(d) detuned by +1T. Due 
to the modulated effective qubit frequency [Eq. (13)], for 7 r /4 the 32 
with detuning —IT is same as the resonant case; for n/2 the 32 with 
detuning —IT is same as the +1T detuned case. The dots in panel 
(d) are the results of the full non-Markovian numerical calculation, 
and the solid curves are based on the Markovian approximation. The 
qubit frequency is too = lOOT. 
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FIG. 14. (Color online) <72 of a single qubit coupled to a semi-infinite 
waveguide with resonant incoming photons. The solid (red) curves 
are based on the Markovian approximation while the (blue) dots re¬ 
sult from the full non-Markovian numerical calculation, (a) koa = 
4l7r/2. Note the breakdown of the Markovian approximation for this 
large value of a. (b) koa = 20n. thus the qubit is at a node of the 
single-photon wave function (a = 10Ao). In the Markovian approx¬ 
imation, the qubit is decoupled from the waveguide and ( 72 (f) = 1 . 
Clearly this is not the case in the full solution—there is both bunching 
and antibunching. The qubit frequency is to 0 = 100F. 


has not been released, and the former can stimulate the emis¬ 
sion of the latter, producing two photons coming out together. 

An additional difference comes from the nodes present in 
the wavefunciton in the semi-infinite case. We find that I = 0 
when kga = 0, 7 r, 27 t, • • •, and hence no bound state is present, 
yielding g? = 1. The qubit, being placed at a node of the 
photonic field, is fully decoupled from the waveguide [49, 90], 

In comparing the koa = 7 r /2 results to those for 7 r/ 4 , it 
is clear that the timescale for features in <72 is larger for the 
smaller value of a. That this should be the case is evident from 
the pole structure: they are symmetric in the n/2 case and 
rotated from that symmetry point for 7 r/ 4 . Thus the lifetime 
for one of the poles in the 7 t /4 case is longer than for the tt/2 
poles, causing the timescale for the structure to be larger. 
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FIG. 15. (Color online) g 2 of two qubits coupled to a semi-infinite 
waveguide with qubit-qubit separation k^L = 7 t/2 . The qubit-mirror 
separation in the first column is koa = 7 t/2 and in the second 7r/4. 
The first row has resonant driving (E/2 = ojo = 100T) and the 
second is detuned by —IT. 


To assess the quality of the Markovian approximation, this 
is one of the cases we have chosen to investigate (for other 
results, see Fig. 6 above). In Fig. 13(d) we compare our an¬ 
alytical Markovian results with the full non-Markovian nu¬ 
merical calculation (blue dots) when the 2LS is very near the 
end of the waveguide, koa = 7r/4. The agreement between 
the two calculations is excellent. However, as koa becomes 
larger, the Markovian approximation breaks down, as demon¬ 
strated in Fig. 14 for koa = 4l7r/2 and 207r corresponding 
to a = lO^Ao and 10Ao, respectively. This seems to hap¬ 
pen when a is larger than a few wavelengths, which for our 
choice of parameters means that a is of order the time of flight 
of a photon during the decay time of the 2LS, a ~ c/2T. 
Fig. 14(b) shows a particularly dramatic example. If the dis¬ 
tance between the 2LS and the mirror is small, placing the 
2LS at a node of the single-particle wavefunction (a = riX (] /2 
for some n) causes it to be completely decoupled from the 
waveguide: there are no incoherently scattered photons, as 
has been discussed theoretically [90] and seen experimentally 
[49], and g 2 (t) = 1 for all t. In contrast, for the large a used 
in Fig. 14(b) so that non-Markovian effects are important, g-> 
shows strong bunching and antibunching. Clearly, the 2LS re¬ 
mains coupled to the waveguide and causes nonlinear bound 
state effects. Though the parameters considered in Fig. 14 fit 
the discussion of non-Markovianity in Ref. [93] in terms of the 
qubit excitation, we leave the problem of a quantitative char¬ 
acterization of the non-Markovianity in this system for further 
study. 

Results for a two-qubit case are shown in Fig. 15. We again 
use koL = 7t/2, giving rise to a symmetric pole distribution 
in an infinite waveguide, and focus on the effect of the mir¬ 
ror. As expected, the oscillation when koa = 7t/ 4 lasts longer 
than that with n/2 due to the existence of the sub-radiant pole. 
This result is consistent with the finding from the calculation 
of power spectra (Figs. 10 and 11). As N increases, we find 
that the behavior of g 2 can be explained in much the same way 
as in the infinite waveguide situation by examining the pole 
distribution. 

We show one representative example of g 2 for 10 qubits and 


resonant driving in Fig. 12(b). While qualitatively similar to 
the result for an infinite waveguide (cf. Fig. 5), g 2 here shows 
a more complex interference pattern and stronger modulation, 
as for the power spectrum. 


VII. CONCLUSION 

In this work we have surveyed a wide variety of multi-qubit 
waveguide-QED structures, focusing on their two-photon non- 
linearities as manifested in the power spectrum [resonance 
fluorescence, ,S'(u;)] and photon correlation function [second- 
order coherence, g 2 (t)]. It is clear that in the multi-qubit case 
(we studied from one to ten qubits), these two functions show 
a great deal of structure caused by the interference of the par¬ 
tial waves scattering from different combinations of qubits. 
Given that oscillations are ubiquitous in g 2 (t) here, the ini¬ 
tial correlation g 2 (t = 0) certainly cannot be used as an in¬ 
dication of whether the system generally causes bunching or 
anti-bunching of photons. 

The structure in g 2 (t) and S(u>) generally becomes sharper 
as the number qubits, N, increases—an effect particularly no¬ 
ticeable in the resonance fluorescence, see Figs. 2 and 4. This 
is natural as the interference effects become more complicated 
and the photonic band gap builds up. In g 2 (t) the deviations 
from semi-classics (g 2 = 1) persist for a much longer time 
than one might initially expect, and this time increases upon 
increasing the number of qubits. For N = 10 the decay of 
correlation in time is very slow indeed (see Fig. 6). 

Many of the features and trends in our results can be roughly 
explained by referring to the poles of the transmission ampli¬ 
tude. These poles (see Fig. 3 for an example) also appear in 
the single particle Green function used in calculating the cor¬ 
relation or “bound state” effects. We have seen that the most 
sub-radiant pole is especially important. The ubiquitous os¬ 
cillations seen come from beating between the frequency of 
the incoming photons (driving frequency) and the real part of 
the most sub-radiant pole. Other oscillations no doubt come 
from beating among the different poles and between them and 
the driving frequency. The long decay time, seen especially 
for large N, comes from the small decay rate of the most sub¬ 
radiant pole; we saw that this scale also appears as the time 
delay. 

Some notable features in our results include: The total 
power spectrum is symmetric about the driving frequency, but 
note that the spectrum of only the transmitted or reflected pho¬ 
tons (in the infinite waveguide case) are not. We have seen that 
there is often either bunching or anti-bunching in both trans¬ 
mission and reflection—because the photons can spend a sig¬ 
nificant amount of time traveling among the different qubits, it 
is not the case that if one is bunched the other should be anti¬ 
bunched. It is unfortunate that there are very few trends as the 
number of qubits increases. One exception is the interesting 
case in which there is strong anti-bunching in transmission and 
bunching in reflection that lasts for a long time (Fig. 7); this 
is enhanced as N increases due to the increasingly sub-radiant 
pole produced by the multiple interference. 

The infinite and semi-infinite waveguide cases show a num¬ 
ber of differences. Perhaps the most important is that a single 
2LS can cause bunching of two photons in the semi-infinite 
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case (Fig. 13) while in the infinite waveguide case there must 
be complete anti-bunching [( 72 (0) = 0]. The mirror in the 
semi-infinite case acts to effectively double the number of 
qubits, and so there is more sharp structure in the presence 
of a mirror for the same number of qubits. 

The effects of loss and dephasing have been entirely ne¬ 
glected in the present study; what effect would they have? Let 
r' denote the rate of decay of one of the qubits to modes other 
than the waveguide. Then, one expects that any structure on 
a timescale larger than (T') -1 will be smoothed out. In par¬ 
ticular, phenomena related to the most sub-radiant pole will 
disappear first, when T' > T. Pure dephasing causes a similar 
smoothing of interference effects without, of course, relaxing 
the excited state population. In addition to smoothing, dephas¬ 
ing can cause the power spectrum to be asymmetric about the 
input frequency [90], One may think, then, that most of the 
structure in our calculated curves would disappear. However, 
there has been tremendous experimental progress recently to¬ 
ward making systems whose loss rate is very low and whose 
dephasing is even smaller. Purcell factors, defined by l /V, 
greater than ten have been demonstrated in more than one ex¬ 
perimental platform: for superconducting qubits coupled to 
a microwave transmission line [43, 45, 46, 49], for instance, 
as well as quantum dots coupled to plasmonic nanostructures 
[68] and to a photonic waveguide [69, 70], Given the rapid 
pace at which the experimental systems are advancing, we 
think considering a system in which P ~ 100 is reasonable; 
for such a Purcell factor, the large majority of the effects pre¬ 
sented here will survive. 

We have considered the effect of weak disorder on our re¬ 
sults, using ±5% as a typical experimental variation in, for in¬ 
stance, the qubit frequency among nominally identically made 
qubits. Such a weak disorder does little to change the results. 
Increasing the disorder would, of course, change the interfer¬ 
ence effects such that the average or typical results would show 
much less structure. However, a large variation in the qubit and 
coupling properties seems unreasonable for the experimental 
systems studied, and in addition, because the system is ID, 
localization of the wavefunctions appears immediately, pro¬ 
ducing very sharp resonances. 

The large majority of our results were obtained in the 
Markovian approximation, as this is the case relevant to 
most current experiments. We have compared to a full 
non-Markovian calculation in a few cases [for example, see 
Figs. 6(d), 13(d), and 14]. If the spacing of the qubits is small, 
then the agreement between the two is very good. However, 
for large separation, we see a substantial difference between 
the two results. For instance, for one qubit coupling to a semi¬ 
infinite waveguide, placing the qubit at a node of the single 
particle wavefunction does not lead to a decoupling of the 
qubit (as it would in the Markovian regime). The criterion 
for when these effects set in is that the separation L should 
be larger than the distance a photon can travel in the decay 
time of a qubit: if L > c/F, a non-Markovian calculation is 
required. Experimentally the non-Markovian regime may be 
reached by connecting superconducting circuits using a coax¬ 
ial cable [1 10], 

We end by simply noting the richness of the correlation phe¬ 
nomena in these waveguide-QED systems and the bright fu¬ 


ture for further investigations because of rapid experimental 
progress. The methods that we use here to study small ensem¬ 
bles of two-level systems could be extended in a straightfor¬ 
ward way to more complex systems, such as three- or iV-level 
systems and photons of different frequency. 
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Appendix A: Matrix elements 


In this appendix we follow the procedure and notation of 
our previous work Refs. [24, 80] to calculate the matrix ele¬ 
ments needed for the (inelastic) power spectra. We focus on 
the calculation for the infinite waveguide cases, and leave the 
semi-infinite cases for the next appendix. Readers interested 
in details of constructing the two-photon wavefunction Ifa) 
should refer to these references. 

As discussed in the main text, the goal is to compute the 
matrix element ct /(</>i(fc)|a ct (:r)|'!/' 2 )- Starting from the single¬ 
particle eigenstate 


\M k ))a 



x)a f R (:r) + </>? (k , x)a[(x)) 


N 




i=l 


|o) 


(Al) 


with the incident photon of momentum k propagating in the 
a-direction, one can construct the two-photon direct-product 
plane-wave state 


|<Mfcl))ai ® \4>l(k2))a 2 ' S j (A2) 

and the associated two-particle identity operator I 2 = 
J2a u a 2 J dk 1 dk 2 \(l) 2 {k 1 ,k 2 )}a 1 ,a 2 {h( k i, k 2 )\- By insert¬ 
ing I 2 into Eq. (4), one observes that the entire two- 
photon wavefunction |^> 2 ) can be expressed using |<j> 2 ), so 
the problem is reduced to calculating the matrix element 
a 1 ( 4 >i{k)\a a (x)\(j) 2 )', the double-momentum integral men¬ 
tioned in the main text is introduced during this insertion of 

By employing the definition of l^i), one observes that 


102(^1; k 2^)ai,a 2 y^2 ^ 


aR(x)\4> 2 {ki, k 2 )} Oil ,0.2 


and hence 


t'R 1 ( k l, X ) 

\M k 2)) 


V2 


Q 2+(2 ^ 1 ) 

(A3) 
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R (0i (k) | a R (x) | fa (fei , k 2 )) 
R(^l(fe)|0R(x)|^2(fcl,fc 2 )) 

R (^l(fc)|o R (x)|^2(fcl,fc 2 )) 

R (^i(fc)|a R (x)|(/> 2 (fci, fe 2 )) 

L(01 (fc) Ia R (x) I <j>2 (fa, k 2 )) 

L{4 , i(k)\a R (x)\4>2(k 1 ,k2)) 

^(<j>i(k)\a R (x)\<f> 2 (ki,k 2 )) 

h(^i{k)\a R (x)\(j) 2 (ki,k 2 )) 


R ' R= 7^ 
UR = 72 
R,L = 7^ 

l,l = 0 
R,R = 0 


L ’ R 72 

RX = 71 

i 

L,L ~ 7! 


^4> R (k 1 ,x)6(k - k 2 ) + (f> R (k 2 , x)S{k - k x ) 
4> R (ki,x)5(k - k 2 ) 

4> R (k 2 ,x)S(k - ki) 

4> R (k 2 ,x)S(k - fci) 

4> R (k 1 ,x)S(k - fc 2 ) 

fi^(fci,a:)(J(fe - k 2 ) + </>R(fc 2 ,a:)(5(fc - £q) 


(A4a) 

(A4b) 

(A4c) 

(A4d) 

(A4e) 

(A4f) 

(A4g) 

(A4h) 


One can compute a ( 0 i(^)l a L(a;)| ( / ) 2 )ai,a, in a similar way. Before giving the final result, we find that defining the following 
four functions (RRi , Ilk,, lAi, , LLi) is useful: 


aPi{k,x) = ^2 


dk[dk' 2 a 


{^i(k)\a^(x )\^ 2 (k[ 

5 ^ 2 )) 011 ,0-2 (MK, k' 2 )\didi) 

E — (k[ + k ' 2 ) + is 


a,(3 = R, L. 


( A5) 


Collecting all the pieces together, the target matrix element 
is given by 

a {<t>i{k)\ap(x)\ip 2 (ki,k 2 ))«R = a {<t>i(k)\ap(x)\<j> 2 (ki,k 2 )) R R 

N 

- "22 a/3i{k,x) (G~ 1 ) l j (djd J \h{ki 1 k 2 )} RR . (A 6 ) 
i,j =1 


From Eq. (A4) and (A 6 ), it is clear that any term multiplying 
the first term above will carry a Dirac delta function, and so 
it will not contribute to the inelastic power spectrum. Only 
the product appearing in the second term can contribute. In 
evaluating that contribution, we find that the following triple 
integrals ( T RR , T LR , T RL , T LL ) simplify the final result: 




aR* ( k , xq) aRj(k, xq + f), 


,0 

T^ l (E,co)= dt e lult 

7—oo 


Jdk aL* (, k, Xo)aLj (fc, xq — t ) 


(A7) 


for a = R,L. We note two things here: First, we call the 
objects triple integrals because there are three momen¬ 
tum integrals to be done; the time integration can be evaluated 
trivially. Second, although formally the function a(3i carries 
an xq dependence, requiring the position xq of the detector to 
be away from the qubit array (xq 0 for transmission and 
io < 0 for reflection), as discussed in the main text, will re¬ 
move this dependence, as expected. 

Finally, substituting Eq. (A 6 ) into Eq. (5) and removing all 
terms proportional to delta functions after integration (they 
represent coherent scattering as discussed in the main text), 
we obtain the incoherent power spectrum 


Sp(u>) = 2 ^ Re [RR( 02 (fct,A: 2 )|d i d i )(G- 1 )* J . (T R0 + T L 2 j k 




(A 8 ) 


where f3 - /f is for the transmission fluorescence and /3 = L for the reflection fluorescence. In the derivation we use the property 
that the transpose of the matrix G is itself. 

To wrap up, we emphasize two things: (a) k\ = k 2 = /s/2 is used, and (b) the Markovian approximation can be introduced 
at the stage of the evaluation of the matrices G and T. For a single-qubit coupled to an infinite waveguide, the double and triple 
integrations can be done exactly, and the result is 

r 4 1 

Sr(w) = Sl(w) = 4^ l(E - w 0 - w) 2 + r 2 /4] [(E/2 - w 0 ) 2 + T 2 /4] [(w - cu 0 ) 2 + T 2 /4] ’ (A9) 

which gives the N = 1 plot in Fig. 2. For N > 2 the results are more complicated, and we do not reproduce them here. 
In particular, for N > 5 it is well-known that there is no explicit formula for the roots of a degree N polynomial, so we use 
Mathematica to symbolically keep track of all the poles and to calculate the measurables presented in this paper. 
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Appendix B: Modifications for the Semi-Infinite Waveguide 


As mentioned in the main text the calculation of the semi-infinite waveguide cases is actually simpler because all the incident 
light can only propagate to the right, be reflected, and then collected at one end. As the first consequence, in contrast to Eq. (Al), 
the single-particle eigenstate \cf>i) no longer carries a directional index a: 


\Mk)) = 


N 

dx £)c4(ir) + (f)L(k, x)a[(x)^j + ^ ei(k)d\ 


| 0 >. 


(Bl) 


Similar to the infinite waveguide cases, to solve for the photon and qubit wavefunctions, ^(k, x), ^(k, x), and (ei(fc)}, the 
following ansatz is used: 


r.ikx 


N -1 


<M fe > x ) = - x) + ^ ti(k)9(x - Xi)6{x i+ 1 - x) + t N (k)0(x - x N )J , 

g—ikx / ^ 

4>h{k, x) = I r(k)9(x 1 - x) + ^ n(k)9(x - Xi)9(x i+1 - x) +r N (k)9(x - x N ) 


(B2a) 

(B2b) 


where r(k) is the single-photon reflection amplitude and 
|r(fc)| 2 = 1 could serve as a check of the calculation. Since 
the number of equations given by the Schrodinger equation 
does not match the number of unknowns, we need one more 
equation to close the set, which is the boundary condition at 
x = 0. In this work, we use hard-wall boundary conditions, 

<fa.(k,0) + 4>h{k, 0) = 0 => t N (k) + r N (k) = 0, (B3) 

where x = 0 is the mirror position. Following the standard 
procedure (see e.g. Ref. [99]) one is then able to compute the 
wavefunctions. 

For the two-photon wavefunction the generalization is 
straightforward. The only change, also stated in the main 
text, is that the single- and double-particle identity opera¬ 
tors, I] and 1-2, carry no directional indexes. In contrast 
to Eq. (A4), in this case the only matrix element needed is 
(</>i(fc)|aL(aj)|<^ 2 (fti, fo))- Therefore, actually one only needs 
to calculate the T RL matrix for the power spectra. This again 
shows the simplicity of the semi-infinite waveguide cases. 

For a single qubit coupled to a semi-infinite waveguide, un¬ 
like the infinite case, however, the Markovian approximation 
is necessary if one wants to analytically evaluate the integrals, 
and the result is the same as Eq. (A9) (up to an irrelevant 
prefactor gone after normalization), except that the qubit fre¬ 
quency and decay rate are replaced by Eq. (13). This result 
agrees with Ref. [90] in the weak driving limit, as we discuss 


next in Appendix C. 

Appendix C: Heisenberg-Langevin Equations: 

One & Two Qubits 

The Heisenberg-Langevin (H-L) equation [102], used in 
Ref. [90], is useful for studying a system under arbitrarily 
large coherent driving. The Langevin (“noise”) term enters 
the equation by integrating out the field degree of freedom and 
carries the information of the driving strength via the initial 
state. In this sense, the H-L equation is formally equivalent 
to input-output theory and gives the same set of differential 
equations [102], 

Specifically, to calculate the measurables using the H-L 
equation, the strategy is to write down the equations of motion 
for both atomic and photonic operators, and then to integrate 
out the photonic degree of freedom, leaving a set of first-order 
differential equations for qubits only. Then one solves the dy¬ 
namics of the qubits so obtains the steady state of the system. 
Substituting the solution back into the formal solution of the 
field gives the photon dynamics. The calculation is explained 
in detail in Ref. [90], so we just list the key changes to arrive at 
the results. Here we focus on one- and two-qubits coupled to 
an infinite waveguide as examples; the generalization to mul¬ 
tiple qubits is straightforward [26, 34] but cumbersome. 

The first step is to get the formal solutions of the photonic 
operators in real space from the Hamiltonian Eq. (1): 


N 

ur(x, t) = clr(x — f, 0) — iV a i~(t — x + Xi)Q(xi < x < t + xi), 

i=t 
N 

a^(x, t ) = a^x +1, 0) — iV cr.;_(i + x — Xi)Q(—t + Xi < x < x{), (Cl) 

*=t 


where the generalized step function is defined by and we have chosen the initial time to be t = 0. Note that 

{ 1, a < x < b 

1/2, x = aorx = b (C2) 

0, otherwise, 









during the derivation it is safe to define the Fourier transform 
of (Ir/l(:e) because of the RWA. 


straightforward, and the result is given by (cf. Ref. [90]) 
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1. One qubit 

Next, one substitutes these formal solutions into the qubits’ 
equations of motion. For a single qubit the calculation is 


(T J ■ 


7 /,A /— Y/2 + i5hj 0 iCl \ f Sl \ 

_ [ a* I = [ o -r /2 - iSu -ifi ■ sj + 

dt \s 2 J V in / 2 ~ iSl / 2 ~ T J \ S 2 / 


(C3) 


where si(f) = (cr_(t)) exp (ikt), s 2 (t) = (<j + (t)<7-(f)), Sco = k — wo, ( 0(t )) represents the expectation value of the operator 
O at time t, and f 1 = \/2TA (note that we use k = E/2 and A to represent the driving frequency ui p and driving strength E of 
Ref. [90]). We emphasize that the equations above can be cast into a compact form: dt S = D S + F with matrix D and vector 
F being constant in time. This /9-matrix is of great importance because it allows one to calculate all higher-order correlation 
functions using the quantum regression theorem (making the Markovian approximation) [ 102], 

Assuming that the qubit is initially in the ground state, one finds si(0) = S2(0) = 0, and so the steady state given by Eq. (C3) 
is 


i r(r /2 + iSuj)n 
Sl(oo) “ 2 r(r 2 /4 + Scj 2 ) + rn 2 /2’ 

( i - 1 

S2looj 2 r(r 2 /4 + Suj 2 ) + m 2 /2 ’ 


(C4) 


Comparing this solution to Eqs. (30) and (31) of Ref. [90], one can see that in general they possess the same form, but there the 
decay rates and the qubit frequency are “renormalized” to the values in Eq. (13) because of the mirror. 

Given the qubit steady states, one can compute the transmission and reflection amplitudes for weak coherent driving. For the 
transmitted part, we take the observation time o > 0 and use Eq. (Cl) to look at the right-going photons: 


(■ a R (x 0 ,T )) = (a R (xo - T, 0)) - iV(&-{T — x o ))0(O <x 0 <T) 

= t f If T(T/2 + iS U )2y/T/2 \ lHxo . T) 
^ 2 V 2 r(r 2 /4 + 5w 2 ) + r 2 A 2 y 


(C5) 


The transmission amplitude for weak coherent driving can be 
defined as 


> _ (&r(T0i T)) A— >0 i 

y ’ ~ Ae ik ^o-T) * 


k — wq 


k — uiq + iT /2 


(C6) 


Similarly, the reflection amplitude for weak coherent driving 

is (x Q < 0,T > |x 0 |) 


n \ = («l(x 0 ,T)) ~iT/2 

U ' } ~ Ae-^o+T) k-ujo + iT/2' 

and the qubit steady state is 


(Cl) 


Let us now turn to calculating the power spectrum. While 
for a single qubit this calculation is straightforward, we sketch 
how to employ the quantum regression theorem to simplify the 
calculation for more complicated cases. Suppose we define 
s 3 (t') = (c + (T)(t_(T + t'))e ikt ', 34 (f) = (<x+(7>+(T + 
f))e- ik(2T+t '\ and s 5 (f) = (a+(T)a + (T + t')<j-(T + 
t'))e ~ lkT , then the quantum regression theorem states that the 
vectors 


S (f) 


/ s 3 \ /-iCls* 1 (oo)/2\ 

I s 4 I , F= I iflsl (oo)/2 I 


(C9) 


_ Si(oo) A->0 

A k — ujq + iT / 2 


(C8) 


These results agree with previous studies [4, 7, 12] (the factor 
v/27r in e(k) comes from the different definition of input state 
for two formalisms). We emphasize that the definitions (C6)- 
(C8) are valid only for weak coherent driving; for arbitrary 
driving amplitude A, in general | (a R ) | 2 + | (ol) | 2 7 ^ A 2 . 


also satisfy dt> S = D S + F. As a result, the matrix D is a sort 
of the characteristic of the system: it gives a closed set of first- 
order differential equations by properly defining the vectors S 
and F. 

Furthermore, in this compact notation we can separate the 
vector S(f) into two parts: S(f) = 8S(t) + S(oo), where the 
former represents the transit dynamics that decays to zero and 
the latter represents the steady state. One can immediately 
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see that the former instead satisfies a homogeneous first-order 
differential equation: dtS S = D6S. If we define a new vector 
I by 


I(w)= / dt' (CIO) 

Jo 

then, after integration by parts, the vector I is given by 

I(w) = - {D + t(w - k)!)- 1 <SS(0), (Cl 1) 

where I is the identity matrix and <5S(0) = S(0) — S(oo). 
This shows that by computing the vector I one can obtain the 
Fourier-transform of the two-time correlations. 

After some algebra we arrive at the transmission power 
spectrum 


-SrM = 


A ~ + 2A \/^ Re ( is *i(°o)) + £ m ° o )| 


S(u> — k) 


+ 2 ^ R- e ^3(w), 


(C12) 


Note that number conservation is indeed guaranteed. 


(orOr) + (a[a L ) = J du (S R (u) + 5 L (w)) = A 2 . (C14) 


Here comes the key point: If one Taylor-expands the inco¬ 
herent power spectrum in terms of the driving amplitude A, the 
lowest-order term is 0(A 4 ), and it would be the one calculated 
from Eq. (3) multiplying by A 4 for dimensional reasons; that 
is, 


S H - L (u) = S L ~ s {w)A 4 + 0(A 6 ) + • • • . (C15) 


Therefore, the two approaches, the H-L equation and the L-S 
formalism, give a consistent result in the weak driving limit. 
In other words, in the weak driving limit the power spectrum 
is solely contributed by two-photon scattering processes [15, 
76]. Note that we have also confirmed this consistency for the 
semi-infinite, single-qubit case considered by Ref. [90], 


where the first term proportional to the delta function repre¬ 
sents the coherent scattering, and the second term represents 
the incoherent scattering. Note that knowing (the com¬ 

ponent of the vector I corresponding to s 3 ) is enough to de¬ 
termine the incoherent power spectra [see Ref. [90] for an ex¬ 
plicit expression for ^(w)]. The reflection power spectrum is 
simply given by 

-SlM = |si(oo)| 2 A(cl; - k) + ~ R.e/ 3 (w). (C13) 

Z Z7T 


2. Two distant qubits 

For two qubits the Markovian approximation is needed for 
an analytical solution, which is also introduced in Ref. [90], 
Based on the single-qubit case, here we list the matrix D de¬ 
rived using the same procedure: 


(iSoj — | 

0 

T p ik 0 L 

2 e 

0 

2 m 

0 

0 

0 

0 

0 

Te ik ° L 

0 

0 

0 

0 ^ 

0 

—iSui — ^ 

0 

r —ikoL 

2 e 

-2 m* 

0 

0 

0 

0 

0 

0 

Y e -^oL 

0 

0 

0 

r ikoL 

2 e 

0 

Uo 

1 

3 

0 

0 

2zor 

0 

0 

0 

0 

0 

0 

TeikoL 

0 

0 

0 

_ r — ikoL 

2 c 

0 

—iSoj — 2 

0 

—2i0l 

0 

0 

0 

0 

0 

0 

0 

p e —ikoL 

0 

m* 

-m 

0 

0 

-r 

0 

_£ e ik 0 L 

Te-ikoL 

0 

0 

0 

0 

0 

0 

0 

0 

0 

im 

-*or 

0 

-r 

p—ikoL 

2 c 

t' ikoL 

2 e 

0 

0 

0 

0 

0 

0 

0 

0 

-zor 

m* 

0 

_r 

_ £ e -ik 0 L 

-r 

0 

0 

0 

-2z9T 

0 

0 

2zor 

2T cos koL 

m 

0 

0 

-im 

I^„—ikoL 

2 e 

p ikoL 

2 e 

0 

-r 

0 

0 

0 

2 m 

-2 im 

0 

2T cos koL 

-im* 

0 

-m 

0 

0 

0 

0 

0 

2i8uj — r 

0 

2z9T 

0 

2 im* 

0 

0 

0 

m 

0 

zOI* 

0 

0 

0 

0 

0 

—2i8uj — r 

0 

-2im* 

0 

— 2 im 

0 

0 

0 

0 

0 

-m* 

0 

-m 

0 

m* 

0 

iSoj — ^ 

0 

r p—ikoL 

2 e 

0 

2 im* 

0 

0 

0 

0 

m 

0 

0 

im* 

0 

-im 

0 

—iSoj — ^ 

0 

r ikoL 

2 e 

- 2 m 

0 

0 

0 

0 

0 

-m 

0 

-m* 

m 

0 

r —ikoL 

2 e 

0 

wh 

1 

3 

(0 

0 

2 im 

0 

0 

0 

0 

0 


z^n 

0 

0 

-im* 

0 

r ikoL 

2 e 

0 

tH|c< 

1 

3 

T 

— 2 im* 

V 0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

m 

-im* 

m* 

-zOT 

-2T ) 


(C16) 
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where 97 = \JY j2 A e lkl > L ! 2 , For equal-time correlations, the vectors S and F are given by 


/Sl (t) \ 

( \ 

/-m\ 

si(ty 


(a 1 + (t))e~ ikt 


m* 

s 2 {t ) 


(& 2 -(t))e ikt 


-m* 

«a(t)* 


(a 2 + (t))e~ ikt 


m 

s 3 {t) 


(<7i + (f)<7i_(t)) 


0 

s 4 (f) 


(a 2 + (t)a 2 -(t)) 


0 

S5(t) 


(a 1 + {t)a 2 _(t)) 


0 

s 5 {ty 

= 

(tT2+(f)cri_(f)) 

, F = 

0 

se(t) 


(ai-(t)a 2 -(t))e 2ikt 


0 

s 6 (ty 


{a l+ {t)a 2 + (t))e~ 2lkt 


0 

s 7 {t) 


(ai + (t)a 1 -(t)cr 2 -(t))e lkt 


0 

S7(ty 


(cri + (f)cri_(f)cr 2+ (f))e- ifct 


0 

Ss(t) 


(ct 2+ ( t)a 2 -(t )<Ji_ (t))e lkt 


0 



(a 2 + (t)a 2 _{t)a 1 + {t))e~ lkt 


0 

\s 9 (t )/ 

\(ai + {t)ai-(t)(j 2 + (t)a 2 -(tD ) 

\ 0 / 


(Cl 7) 


Following the same procedure, one can construct the two- and multiple-time correlations using the same matrix D. We checked 
that the results agree with those of Ref. [26], as well as those obtained from the L-S formalism presented in the main text. Finally, 
we note that for N qubits in the Markovian regime, the size of the matrix D grows as (4' v — l) 2 , where the factor of 4 arises 
because each qubit can be represented by a xtVtZ or 1 and subtracting one means not to consider the trivial correlation (1 :V ). 


Appendix D: Two-photon Transmission Probability 


In this appendix we consider the transmission and reflection probabilities using the L-S formalism. Specifically, we calculate 
the expectation values (V^IorOrI^) and . First °f all, this calculation serves as a consistency check, and it can 

also be compared with the H-L approach. Secondly, it justifies our argument that it is not necessary to consider the transmission 
<72 when driving the system resonantly. Here we focus on two illustrative examples: a single qubit coupled to both infinite and 
semi-infinite waveguides. 

We start from the transmission probability for the infinite waveguide case. One first inserts the one-particle identity operator 
Ii and employs Eq. (A6), 


{4>2\a f R (x 0 )a R (x 0 )\ip2) = I dk[ \ R {(t>l{k)\a R (x 0 )\il 2 {ki,k 2 ))\ 2 +\l J {(t) 1 {k)\a R {xo)\'^l) 2 {k 1 ,k 2 ))\^ 

pik ix 0 


/H( : 


\Z4n 


t(k\)5{k — k 2 ) + ki 0 k 2 I — RR(k, Xq)G 1 e(ki)e(k 2 ) 


+ / dk 


LR(k,x 0 )G 1 e(fci)e(fc 2 ) 


(Dl) 


One can see that the above expression contains three parts: the plane waves, the bound state, and their interference. The integra¬ 
tion over the plane waves is 


/ 


dk 


gifeixo 

— j^^t{ki)5(k — k 2 ) + ki <+ k 2 
V47T 



<5(0) + 


t{kiYt{k 2 ) 

-7+-<H fc i 



+ k\ 7A k 2 , 


(D2) 


and the interference term is 


- / dk 


/ e ik 1 x 0 

f _ t(ki)S(k — k 2 ) + fci <+ k 2 ) x RR(k,x 0 )G~ 1 e(ki)e(k 2 ) 


— h.c. = 0. 


(D3) 


Finally, the bound state term is 


dk 


\RR(k,x 0 )G 1 e(fci)e(fc 2 )| +1 LR(k,x 0 )G 1 e(fci)e(fc 2 )| 


2 r 3 /r 


[(E - 2oj) 2 + T 2 ] 


2 > 


(D4) 


which could also be derived by integrating Eq. (A9) over to 
and dividing by 27 t. Combining the three pieces together and 


taking k\ = k 2 = E/2, we have 

<^l4« R |^2) = l^(-s/2)i 2 — + — YUM. — 

R ' 7 T l(E - 2w) 2 + r 2 ] 2 

(D5) 
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which is indeed independent of xq as it should. 

The presence of the infinity S(k = 0), proportional to the 
“volume” of the system, in Eq. (D5) seems rather awkward and 
needs interpretation. Our one-photon input state, implicitly 
defined in Eq. (B2), is 



which gives (k\k') = 5{k — k') and (fc|a^aR|fe) = l/2n. 
Using this to construct the two-photon input state, \k \, k 2 ) = 
|fci) ® \k 2 )/y/2, one obtains 


{ki,k 2 \a^a R \ki,k 2 ) = 2 x —— (D7) 

Z7T 

for k\ = k 2 . This is dimensionally correct and tells us that 
the “number of photons” injected into the system is 8(0)/n. 
Therefore, dividing Eq. (D5) by S(0)/n gives the two-photon 
transmission probability 


T 2 


(ip2\aiaR\4’2) 
8(0)/it 


IW2)| 2 + 


4 r 3 1 

[(E ~2co) 2 +r 2 } 2 2ttS(0)' 


(D8) 


enough (A 1 2 <C T), the two-photon correction to the transmis¬ 
sion probability is negligible and we do not have to consider 
the transmission g 2 when driving resonantly (so T = 0). 

One could do the calculation of (il) 2 \a[ai\'^ 2 ) in a similar 
way. The three pieces (plane wave, interference, bound state) 
are 


r(E/2)\ 


2 m 


-AT 3 4 5 6 7 /? 


2T 3 /7 


7T ' [(£ - 2w) 2 + T 2 ] 2 ’ [(E-2u) 2 + T 2 Y 

(D9) 

respectively. We note that in contrast to the transmission, the 
interference term in reflection is not zero and that the last term 
(bound state) is the same as that for transmission. Therefore, 
the two-photon reflection probability is given by 


I 2 ! 


r 2 


<5(0)/7t 


HE/ 2)| 2 


4T 3 1 

[(£-2w) 2 +r 2 ] 2 27n5(0)' (D10) 


One can see that indeed T 2 +R 2 =T + R=1 holds, which 
relies on the precise interference between the plane wave and 
bound state parts. 

This precise interference can also be seen in the semi¬ 
infinite waveguide case. Without transmission, one only needs 
to compute (ip 2 \alai/tl) 2 ) and the three terms are given by (in 
the Markovian regime) 


This dimensionless expression is transparent: the first term is 
the single-photon transmission probability T = \t(k)\ 2 , and 
the second term is the two-photon correction. Compare this 
expression with that calculated from the H-L equation (cf. Ap¬ 
pendix Cl); in the weak driving limit, the results agree upon 
requiring A 2 = 1 /27t< 5(0). Therefore, if the driving is weak 


2 5(0) -16T 3 /7t 2 


\r(E/2)\ 2 ^, - 


16f 3 /7 


(E - 2w) 2 + f 2 (E- 2(D) 2 + f 2 


(Dll) 

As a result, the last two terms cancel exactly, leaving R 2 = 1. 
This completes the consistency check of our theory. 
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